Existence, Stability, and Scattering of Bright Vortices in the Cubic-Quintic Nonlinear 

Schrodinger Equation 
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We revisit the topic of the existence and azimuthal modulational stability of solitary vortices 
(alias vortex solitons) in the two-dimensional (2D) cubic-quintic nonlinear Schrodinger equation. 
We develop a semi-analytical approach, assuming that the vortex soliton is relatively narrow, and 
thus splitting the full 2D equation into radial and azimuthal ID equations. A variational approach 
is used to predict the radial shape of the vortex soliton, using the radial equation, yielding results 
very close to those obtained from numerical solutions. Previously known existence bounds for the 
solitary vortices are recovered by means of this approach. The ID azimuthal equation of motion is 
used to analyze the modulational instability of the vortex solitons. The semi-analytical predictions 

- in particular, that for the critical intrinsic frequency of the vortex soliton at the instability border 

- are compared to systematic 2D simulations. We also compare our findings to those reported in 
earlier works, which featured some discrepancies. We then perform a detailed computational study 
of collisions between stable vortices with different topological charges. Borders between elastic and 
destructive collisions are identified. 



I. INTRODUCTION 



ferred to as its topological charge (alias vorticity), m. An 
example of a stable solitary vortex is depicted in Fig. Q] 



The cubic-quintic nonlinear Schrodinger (CQNLS) equa- 
tion is used to model a variety of physical settings. In 
scaled units, the CQNLS equation takes the following form: 



v 2 * 



I*! 2 * - Ifl 4 * = o, 



(i) 




where ^>{x,y,t) is the complex wave function, V 2 is the 
two-dimensional (2D) Laplacian, and the last two terms 
represent, respectively, the focusing cubic and defocusing 
quintic nonlinearities. The CQNLS equation emerges in 
models of light propagation in diverse optical media, such 
as non-Kcrr crystals [2fJ, chalcogenide glasses 
ganic materia ls J3H , colloids @, Q , dye solutions 
ferroelectrics [17| ■ It has also been predicted that this com- 
plex nonlinearit y ca n be synthesized by means of a cascad- 
ing mechanism [ljj. It should be noticed that, in the op- 
tics models, evolution variable t is not time, but rather the 
propagation distance. 

The competition of the focusing (cubic) and defocusing 
(quintic) nonlinear terms is the key feature of the CQNLS 
model, which allows for the existence of stable multidimen- 
sional structures which would be unstable in the focusing 
cubic nonlinear Schrodinger (NLS) equation. In particu- 
lar, the CQNLS equation supports solitary vortices (alias 
vortex solitons) in 2D and 3D geometries. These are ring- 
shaped structures which carry a rotational angular phase. 
The respective integer winding number in the vortex is re- 
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FIG. 1: (Color online) An example of a vortex soliton with 
charge m = 5 and frequency Q = 0.16, found as a numerical 
solution to the CQNLS equation. The height of the profiles rep- 
resents the value of the wave function. The gray profile displays 
the squared absolute value, while the dark (blue) and light (red) 
meshes represent the real and imaginary parts, respectively. 

It is well known that the cubic NLS equation supports 
"dark" (delocalized) and "bright" (localized) vortices, with 
the self-defocusing and focusing signs of the nonlinearity, 
respectively. Dark vortices of unitary charge are stable, 
while higher-order ones are unstable, splitting into unitary 



eddies [13|. On the other hand, the "bright" modes (vortex 
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solitons) are always azimuthally unstable in the framework 
of the cubic equation or its counterpart with the saturable 
nonlinearity [14| , which leads to the breaking of the vortex 
soliton into a number of fragments depending on the charge 
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[H [24[ . The addition of a periodic potential to the cubic 
NLS equation may stabilize vortex solitons of a different 
type, which are built (in the simplest case) as a chain of four 
peaks with the phase circulation corresponding to integer 
vorticity to [!, [H, HH, H(| H3] • An example of the azimuthal 
breakup of an unstable vortex is displayed in Fig. 
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FIG. 2: (Color online) An example of the evolution of the vortex 
soliton with charge m — 4 in the CQNLS equation, exhibiting 
breakup (left to right) due to the azimuthal modulational in- 
stability. Shown is the intensity (squared absolute value of the 
wave function) at consecutive moments of time, from left to 
right. The numerical method and parameters are presented in 
Sec. El 

In contrast to the cubic NLS, the CQNLS can support 
stable solitary vortices. For the first time, this remarkable 
fact was discovered "empirically" in simulations of the 2D 
CQNLS equation [3l|, and later was investigated in detail 
in a more rigorous form in Refs. [23l. |36|. Moreover, three- 
dimensional solitons with embedded vorticity to = 1 also 
have their stability region in the 3D CQNLS eauation[2r| 
(see also reviews [ljj, [24|)- Stable vortex solitons may find 
their potential applications in the design of all-optical data- 
processing schemes. In that respect, knowing the growth 
rates of unstable modes is also important, because, if the 
rates are small enough, the vortices may be considered as 
practically stable ones, as they will not exhibit an observ- 
able instability over relevant propagation distances. 

One of objectives of this paper is to revisit the topic of 
the azimuthal modulational stability of solitary-vortex so- 
lutions in the 2D CQNLS equation. As mentioned above, 
numerous studies of this problem have already been pub- 
lished [hJ [H H3, IM HI. The aim of the analysis was to 
find a critical value (f2 s t) of the intrinsic frequency of the 
vortex [for its exact definition see Eqs. @ and below], 
above which the vortices with a given value of the topolog- 
ical charge are stable. In Ref. [31(, 2D azimuthally stable 
vortices with to = 1 were shown to exist. It was found that 
the slope of the vortex' profile at the pivotal point peaked 
at a specific value of f2, which was considered as the crit- 
ical frequency, Q. s t(m = 1) w 0.145, while the vortices, 
with all values of to, exist at 17 < fi, 1 ^ = 3/16 = 0.1875 
(at f2 = f2„^ x , the radius of the vortex diverges, i.e., 
the "bright" vortex goes over into a "dark" one). The 
same value, f2„^ x is simultaneously the largest one up to 
which exact soliton solutions exist in the ID version of the 
CQNLS equation (28[. In work [3l|, a variational approach 
(VA) was developed for 2D vortex solitons, yielding results 
similar to those obtained in the numerical form. Full 2D 



simulations of the CQNLS equation reported in Ref. [3l| 
had confirmed that vortices with f2 st < O < S!, 1 ,^ were 
indeed stable, while those with < O < fi 8t were not. 

Works [13, [H, Hl| employed a more rigorous approach. 
They introduced small perturbations around the 2D vor- 
tex soliton and solved the resulting eigenvalue problem nu- 
merically. In this way, the critical frequencies were found 
as n Bt (m = 1) « 0.16 and O st (m = 2) « 0.17. Also, 
in Ref. [29[, using the Gagliardo-Nirenberg and Holder in- 
equalities together with Pohozacv identities, is was shown 
that the eigenvalues generated by the CQNLS equation 
possess an upper cutoff value. 

In Ref. 27J, 2D perturbations were considered too. 



Through an extensive analysis, the problem was trans- 
formed into finding, by means of numerical methods, ze- 
ros of respective Evans functions for a set of ordinary dif- 
ferential equations (ODEs). In this way, the existence of 
azimuthally stable vortices of all integer values of m was 
predicted (in Refs. [36[ and [3l[, the stability regions for 
to > 3 were not identified, as they are extremely small). 
The predictions for the critical frequencies reported in Ref. 
j27| . which are somewhat different from those in Refs. [36| 
and [3l| , are shown in Table HU 

The present work aims to undertake an additional study 
of the stability of vortices in the CQNLS equation for two 
reasons. Firstly, the previous studies demonstrated some 
(relatively small, but not negligible) disagreement in the 
predictions of the critical frequency. Therefore, since we 
will use a different approach, our results may help in com- 
paring and verifying the previous findings. Secondly, al- 
though in Ref. [3l[ 2D simulations were performed for the 
vortices, this was only done for m = 1, and longer simula- 
tions later showed that some vortices which were originally 
concluded to be stable turned out to be eventually unsta- 
ble [36|. Therefore, stability results from additional 2D nu- 
merical simulations are needed for comparison with various 
predictions of the stability. In fact, our numerical results 
will illustrate that the predictions of Ref. (2?| are the most 
accurate ones. In parallel to extensive direct simulations, 
we elaborate a semi-analytical methodology, following the 
lines of previous studies of the azimuthal modulational in- 
stability of vortices in the cubic NLS equation Q. The 
approximation amounts to splitting the full 2D equation 
into effectively one-dimensional radial and azimuthal equa- 
tions. The former equation is used to compute the radial 
shape of the soliton by means of a VA to find an analytic 
approximation to the profile, which is then used as an ini- 
tial condition to a numerical optimization routine to find 
the numerically 'exact' profile. The closeness between the 
numerically computed profiles and the VA approximation 
testify to the good accuracy of the VA. The analysis of 
the azimuthal equation makes it possible to predict the 
threshold of the onset of the splitting instability of the 
vortex soliton. Our results predicting the growth rates of 
individual azimuthal modes of unstable vortices match full 
simulations very well, but the predictions for the critical 
frequency are less accurate, when compared to the Simula- 
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tions. This semi-analytical technique may be quite relevant 
for applications in other 2D models. 

Another objective of the work is to study, by means of 
systematic simulations, collisions between stable vortices. 
Except for few examples reported in Ref. [3lj , this problem 
was not studied before. 

The paper is organized as follows. In Sec. 2 we derive an 
approximate analytical description of the vortex profile by 
first finding an analytical asymptotic approximation to it, 
which is then employed as the ansatz on which the VA is 
based. In Sec. 3 we use the variational ansatz as the initial 
guess, to generate numerically exact vortex profiles by dint 
of a nonlinear optimization routine. We then compare the 
numerical profiles with the variational ansatz to show its 
accuracy. In Sec. 4 we derive an approximate ID equation 
for the dynamics along the azimuthal direction. A linear 
stability analysis is subsequently performed to find stabil- 
ity criteria and growth rates of unstable modes within the 
framework of the azimuthal equation. In Sec. 5 we present 
full 2D simulations of the vortices, and the respective re- 
sults for their stability. These results are compared to our 
predictions, as well as to the predictions produced by pre- 
vious works. In Sec. 6, we use direct simulations to explore 
collisions between stable vortices in detail. In particular, 
these studies allow us to identify a border between quasi- 
elastic and destructive collisions. Finally, in Sec. 7, we 
summarize our findings and formulate concluding remarks. 

II. APPROXIMATE ANALYTICAL PROFILES OF 
STEADY-STATE VORTICES 

A steady-state vortex solution to Eq. ([1]) is looked for as 
*(r t 9,t) = f(r)A(9,t), (2) 

where real function f(r) is a stationary radial profile with 
azimuthal dependence given by 

A(6,t) =e t{me+nt \ (3) 

with topological charge m and frequency tt. Analytical 
solutions being not available for the profile f(r), we begin 
the study by finding an approximate analytic expression for 
f(r) and identifying its existence bounds. As shown below 
in Sec. IIII1 the predicted profile is very close to the true 
solution, therefore it can be used to predict the existence 
and stability regions of the vortex solitons. The developed 
analytical method is quite general and may be applied to 
other physically relevant equations. 

A. The Asymptotic Vortex Profile 

Inserting expression (J2J into the underlying equation (JTJ 
yields the following ODE for the radial profile of the vortex 
with charge m: 

~T ( r f) ( n + ^t) /W + / 3 « - = °- ( 4 ) 

r ar \ ar J \ r A J 



If we assume that the vortex takes the form of a relatively 
narrow ring with a large radius, then, in the region of in- 
terest, variable r may be approximately replaced by a con- 
stant, r c 3> 1, which we take to be the value of r at a point 
where the radial profile attains its maximum. Since we as- 
sume r c to be large, in the lowest approximation we neglect 
the 1 jr c term in the Laplacian. An obvious consequence of 
dropping this term is that the approximate profiles which 
we are going to derive will be less accurate for smaller val- 
ues of r c , see Sec. Mil 

Under the large-radius assumption, Eq. (@| becomes 

0-fr/M + / 3 (r)-/ 5 (r)=O, (5 ) 

n* = n + m 2 /r 2 c . (6) 

For a relatively narrow ring, one can treat Eq. ([5]) as if 
r G (— oo,+oo), in which case the equation has a well- 
known ID analytical solution for the soliton [9l l28j: 



1 + y/l - (16/3)0* cosh ( 2y/Q* (r - r c )j 

As mentioned above, this solution exists for < O* < 
fi^ = 3/16. At point O* = 3/16, the solution degener- 
ates into a constant, f(r) = fa = \/3/4. The family of 
solutions ([7]) parameterized by 0* is depicted in Fig. [3] 
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FIG. 3: (Color online) A set of solutions J7j for D.* = 
0.05,0.1,0.15,0.18,0.1874,0.1874999,0.1875 are shown from 
bottom to top. As fl* increases, the profile flattens out, and 
it degenerates into a constant at Q — Q* . 

To obtain a closed-form analytical approximation, it is 
now necessary to find an expression for the radius of the 
profile, r c , for given frequency Q. To this end, we will use 
the VA, and will then compare the results with numerically 
found profiles. 
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B. The Variational Approach 

To apply the VA, we use the Lagrangian density corre- 
sponding to the CQNLS equation: 

C = l - (*** - tf*tf t ) + |* r | 2 +r- 2 \^ e \ 2 - i|v]/| 4 + i|3f, 

(8) 

where the subscripts denote partial derivatives. Inserting 
here expression @, to be used as a factorized ansatz, yields 



and the other one which is equivalent to Eq. ^ . Eliminat- 
ing r c from these equations, one concludes that f2* is, for 
a chosen value of fi, a solution to the following equation, 



o ^ 3 

n = — i 

2 32 16T 



= G(ST). 



(15) 



This yields the final form of the vortex radial profile, as 
predicted by the VA: 



4ft* 



df 



1 



1 



C ={ Q + ^)^ + [Tr) -iArJ+gAr). (9) 

The integration of density ((9]) gives rise to the full La- 
grangian: 

r°° / ii 

L = 2n C{r)dr = 2tt I il d + C 2 + m 2 C 3 - - C 4 + - C 5 

(10) 

where 

d= J™ f 2 (r)rdr, C 2 = J™ ($L\ r dr, 



1 + ^/l - (16/3)0* cosh ( 2^/^F (r - r™) 



(16) 
(17) 



C 3 
C 5 



f(r)-, C 4 



f(r)rdr. 



f\r)rdr, (11) 



Finally, the transcendental equation (fT5")) . together with 
Eq. (|13[) . was solved numerically by means of a simple bi- 
section method with a tolerance of 10~ 15 . 
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Next, we use the asymptotic approximate profile (J7JI as 
an ansatz for radial profile, treating i7* and r c in Eq. (JT)) 
as variational parameters. Inserting expression (0 into the 
integral terms in the Lagrangian, and again making use of 
the large-radius approximation yields 
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(12) 



— VfF-3^T^— -n* 
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T = arctanh 



(13) 



The respective static Euler-Lagrangc equations, dL/dr c 
and BL/dVt* = 0, yield, respectively, equation 



n — 



/3H* 



-1/2 



16 8T 



(14) 



FIG. 4: (Color online) Top: The vortex' radius, as predicted 
by the variational approximation, r™, versus = G(Q*), and 
Q*, for m = 5. Bottom: The values of f2 = G(fi*) vs. fl* , see 
Eq. @. G(fi*) was evaluated over interval fi* £ [0.01,0.1875] 
with step AO* = 0.00001. 

In Fig. 2] we show versus f2 and f2* for m = 5, along 
with relationship (|15p . We see that, with the increase of 
f2, the radius of the vortex starts from very large values 
(infinite at fi = 0), decreases until it reaches a minimum 
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value, and then increases rapidly (it becomes infinite once 
again at fl = fl™ x ). We also see that the relationship 
between fl and fl* does not depend on charge m, and that 
fl — > 3/16 as fl* — > 0.1875 (the apparent gap near the right 
edge of the plot is due to the sensitivity of the relationship, 
which we discuss below). 



C. Existence Bounds for the Two-Dimensional 
Vortex Profiles 

From numerical investigations performed in Rcfs. (3l| 
and [ll|, the existence border for the two-dimensional 
CQNLS equation vortex solution was found to be 
(^max) num ~ 0.180, while, as said above, the analyti- 
cal limit, which is identical for ID and 2D equations, is 
= ^iL = 3/16 = 0.1875. The discrepancy may be 
explained by the sensitivity of relationship (|15p between 
fl and fl*, as produced by the variational approximation. 
From Eqs. $T5§ and (H3J) we see that, as fl* 3/16, 
T -> co, and so fl = G(fl*) -> fl*/2 + 3/32 = 3/16, in 
accordance with the exact result. However, as seen in the 
right panel of Fig. @J the relationship between fl and fl* 
seems to be extremely sensitive near this limit (see the in- 
set in the right panel of Fig. [4]). This sensitivity can also be 
observed in Table HI where some values of fl* and the corre- 
sponding G(fl*) = fl are given. It is clearly seen that one 
quickly approaches the limit of machine double precision 
for fl* as fl -> 0.1875. 



fl* 


n = G(fl*) 


0.1874 


0.1664... 


0.187499 


0.1736... 


0.187499999 


0.1783... 


0.187499999999 


0.1806... 


0.1874999999999999 


0.1823... 



TABLE I: The evaluation of fl = G(fl*) near = 0.1875 

for double precision arithmetic. Four significant digits in G(fl*) 
are given. The extreme sensitivity of relationship fl — G(fl*) 
is observed. To obtain more precise fl* values for fl closer to 
flSL = 0.1875, the use of higher-precision arithmetic is re- 
quired (results not shown here). 

This effect can also be understood from the logarithmic 
divergence of T(Q*) close to the limit point, see Eq. (|13l) . 
Therefore, it is not surprising then that numerical esti- 
mates of fl, 2 ,^, obtained by means of a shooting method 
to look for profiles at different vales of fl, gave lower es- 
timates of the existence bound. In Sec. IIII1 we confirm, 
using accurate numerical methods to find vortex profiles 
for A > 0.180, that the actual existence bound flj ^L i s 
greater than the numerical estimates of Refs. [3l[ and jllj. 

The analytical approximation (| of the vortex profile 
is used in the following sections as an input to a numeri- 
cal optimization routine which solves the ODE ((4]) to find 



numerically exact radial profiles. 



III. NUMERICALLY EXACT STEADY-STATE 
VORTEX PROFILES AND COMPARISON TO THE 
VARIATIONAL APPROXIMATION 

To find numerical solutions to Eq. (|4]), we used a modi- 
fied Gauss-Newton optimization routine, with a tolerance 
10~ 7 [8j. In Fig. [5j some numerically found radial pro- 
files are displayed for charges m = 1, m = 2, and m = 3 
and fl > 0.18, including fl as close to = 0.187 5 as 

fl = 0.186. To our knowledge, the profiles above fl = 0.181 
have not been shown in any previous study. It is seen that 
the profiles flatten out as fl grows, pushing the profile far- 
ther from r = 0. At the same time, the ring of the vortex 
becomes wider without a change in its height. It is worth 
mentioning that, due to the high sensitivity of the numer- 
ics mentioned in the previous sections, computing profiles 
past fl = 0.18 is a daunting task requiring high-precision 
arithmetic. Implementing high-precision arithmetics in a 
Gauss-Newton optimization routine (or any fixed-point it- 
eration method) would be quite involved and very time 
consuming. Nonetheless, we have found that, using high- 
precision arithmetics (up to 300 decimal places for frequen- 
cies as close to 0.187 5 as 0.187) in the calculation of fl* 
in the variational equation (|15[) yields approximate analyt- 
ical solutions that the Gauss-Newton subroutine is able to 
process into numerically accurate solutions for values of fl 
very close to fl 2 ^ = 0.1875, as shown in Fig. [5] The re- 
sults clearly demonstrate the usefulness of the VA: without 
using this approximation for generating the initial profiles 
to be fed into the numerical solver, a prohibitively com- 
plex high-precision fixed-point algorithm would be needed 
to obtain meaningful profiles past fl = 0.18. 

The accuracy of the variational ansatz per se can be 
tested by comparing it to the numerical solutions computed 
by means of the Gauss- Newton routine. In Fig. [6] we com- 
pare the VA prediction for the vortex-ring's radius, r c , to 
the numerically found values for m = 5. We observe that 
the radii provided by both methods are virtually identical, 
allowing one to use the VA-predicted radius in applications, 
that may be useful in experimental situations. 

To get an even better idea of how close the variational 
approach is to the Gauss-Newton profile overall, in Fig. [7] 
we compare the relative sum of squared deviations of the 
variational profiles from their the numerically found coun- 
terparts, for different values of m and different values of fl. 
We also plot a series of profiles produced by the VA and by 
the Gauss- Newton method for fl = 0.15. We observe that, 
as m increases, the mismatch between the variational and 
numerical profiles decreases. This is understandable, as we 
have neglected the l/r c term in the Laplacian, while for- 
mulating the variational ansatz. Therefore, since vortices 
with smaller m have smaller inner radii, the discrepancy 
between the variational ansatz and the numerically exact 
solution are expected for lower values of m. The total dis- 
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FIG. 5: Steady-state vortex radial profiles computed by means of the Gauss-Newton routine for m — 1 (top left), m = 2 (top 
right), and m — 3 (the bottom row of panels) for the indicated values of Q. In all cases, the variational ansatz was used as the 
input (not shown here) with the value of ST computed using high-precision arithmetic (see text for details). The radial direction 
was discretized with a grid spacing of Ar = 0.5. The stopping tolerance of 1CF 7 was used in the Gauss-Newton routine. As O 
approaches the limit value, 3/16, the profiles flatten out due to the increase of the vortex' radius. 



crepancy is observed to be quite low overall, showing that 
the variational profile provides for a very accurate rendi- 
tion of the true solution, especially for larger values of m. 
Since our VA provides an accurate radial profile of the vor- 
tex solitons, we can use it to derive fully analytic azimuthal 
modulational stability predictions, which is done below. 



IV. THE AZIMUTHAL MODULATIONAL 
STABILITY: ANALYTICAL RESULTS 

With the profiles of the steady-state vortices available, 
we proceed to the study of their azimuthal modulational 
stability. To this end, we apply the methodology utilized 
in previous works, dealing with circular gap solitons in the 
model of the circular Bragg grating Q, and in the cubic 
NLS equation Q: we will derive an azimuthal equation 
of motion by assuming a frozen radial profile, and then 
perform a perturbation analysis to analyze the stability of 
azimuthal perturbation modes. 



In Ref . [|| , numerically computed eigenmodes of pertur- 
bations around the solitary vortices featured a slight cou- 
pling in the radial and azimuthal directions, hence assum- 
ing the compete separability of the radial and azimuthal 
directions for the evolution of the perturbations may lead 
to discrepancies between the predictions for the stability 
and numerical results. Nonetheless, the coupling between 
the radial and azimuthal directions is attenuated with the 
increase of m, hence, the predictions made under the as- 
sumption of the separability should be more accurate for 
higher values of m. The trend to the improvement in the 
accuracy of the analytical approximation with the increase 
of m was observed in the numerical results presented in 
Sec. HIl 



A. The Azimuthal Equation of Motion 

To derive the azimuthal equation of motion, we first in- 
sert the separable ansatz ([2]) into the Lagrangian density 
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FIG. 6: The vortex-ring radius r c versus Q as predicted by the 
variational approximation (r™, white circles), and produced by 
the Gauss- Newton profile (rf n , black squares) for m = 5. The 
numerical radius was computed as the radial center of mass of 
the vortex ring. 



to obtain 



£ = zf(r) (AA* t - A* A t ) 



(18) 



l/ 2 (r)|V-l/ 4 (r) |A| 4 + l/ 6 (r) \A\< 



Using our steady-state radial profiles, we can perform the 
integration of the Lagrangian density over dr, thus arriving 
at an effectively one-dimensional (in the direction of 9) 
Lagrangian that may be used to derive the equation of 
motion for A(8, t): 



£ 1D = -CMA* t 



£id do, 



A*A t ) + C 2 \A\' 



(19) 



(20) 



+C 3 \A e 



si 



Ci\A\ 4 



S2 



C 5 \A\ e , 



and Cj, j = 1, 5, are the radial integrals defined as per 
Eq. (|TTj) . Evaluating the variational derivative of the action 
functional [22j . and applying a linear transformation, 

A -> A exp^-idt/d), dt -> Cat/Cx, (21) 

yields the following evolution equation for A(8,t): 

iA t = -A gg - (C4/C3) I A\ 2 A + (Os/Cs) I A| 4 A (22) 

This equation provides a description of how the radially- 
frozen, azimuthally time-dependent solution will evolve in 
the CQNLS model. 
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FIG. 7: Top: The sum of squared errors between the varia- 
tionally predicted and numerically found radial profiles of the 
vortex soliton. Bottom: Examples of the variational (dashed) 
and numerical (solid) profiles for Q — 0.15. Shown from left 
to right are the profiles for m = 1,...,5. It is seen that the 
error decreases as m increases. The series of error values for 
Q = 0.05 and 0.1 terminate at m — 8 because at that point the 
variational profiles are already within the prescribed tolerance 
of 10" 6 of the Gauss-Newton routine. 



B. The Stability Analysis 

To study the azimuthal modulational stability of vortex- 
soliton solutions to the CQNLS equation, we performed a 
perturbation analysis in the framework of the azimuthal 
equation of motion (|22p . with the objective to compute 
the growth rates of small perturbations. We begin with an 
azimuthal "plane-wave" solution perturbed by a complex 
time-dependent perturbation: 

A{6,t) = [l + u(9,t) + iv{9,t)} e l ( me + n ' 4 ), (23) 



with \u\, \v\ <C 1. Inserting this into Eq. (f2"2"j). performing 
the linearization with respect to the perturbations, and 
separating the result into real and imaginary parts, we ob- 
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tain a system of coupled equations (a time shift was made 
here, t t + 6/ (2m)): 



Ui 



-vee, 



v t = u es + I — 



(24) 



In order to study the azimuthal modulational stability of 
the solitary vortices, we expand u and v into Fourier series 
over azimuthal harmonics with integer wavenumbers K: 



u(K,t) 



v{K,t) 



2jv 



u(6,t) e lKti dO 



2jv 



v(0,t) e lKe d9 



(25) 



Applying these transforms to Eq. (f2"4")l yields two coupled 
equations for the amplitudes of u and v of each mode. In 
a matrix form, the equations are 



d_ 

dt 







u 




V 







2C 4 - 4C 5 
C~ 3 



-K- 



K 2 




(26) 



The growth rates for each azimuthal wavenumber K are 
simply eigenvalues of Eq. (|26|) . Taking into account the 
underlying rescaling (|21[) . they are: 



(27) 



where -Kcrit is the critical value of K, above which all the 
modes are stable: 



Kc.rit = 



2(C 4 -2C 5 



C, 



(28) 



For the modulationally unstable vortices, it is useful to 
know the maximum growth rate and the mode that exhibits 
it. This is because, in an experiment, even an unstable 
vortex may be practically "stable enough" if the maximum 
perturbation growth rate is small enough. This informa- 
tion can be extracted from expression (|27[) by equating 
its derivative to zero and solving for K , which reveals the 
fastest growing perturbation mode (and subsequently, the 
prediction of the number of fragments that the unstable 
vortices will break up into): 



K„ 



C*4 — 2C*5 



Co 



V2 



(29) 



We can then insert this value into Eq. ([28]) . which yields 
the maximum growth rate, 



A, 



Ci ~2Ch crit 



(30) 



In fact, because K must be integer, the actual fastest grow- 
ing eigenmode may correspond to the integer K closest to 
value (|2"5|) . Accordingly, the actual largest growth rate may 
be somewhat smaller than the one given by Eq. (J3UJ) . 

For the vortex to be azimuthally stable against all modes, 
one needs either K CI n < 1 (then, there is no integer value 
K < -Kcrit) or -Kcrit being purely imaginary. According to 
Eq. (|2"5)l . these two stability criteria amount to the follow- 
ing inequalities: 

imaginary K cl - it : C 4 - 2C 5 < 0; (31) 

|tf cri t| < 1 : Ci - 2C 5 - C 3 /2 < 0. (32) 

Since C3 is positive, the former condition is stricter than 
the latter one. Although the latter condition is sufficient 
for the azimuthal stability, we also keep the former one 
(according to the above derivation, it is relevant to the in- 
finite line) because it leads to an estimation for the critical 
frequency which is independent of charge m (see below). 

For our predictions, we compute the Ci constants using 
numerical integrals of the numerically exact radial profiles. 
However, since the variational analytic profiles are so close 
to the numerically exact profiles, one can use expressions 
(fl"2)) to find approximate analytic values of Eqs. (f3"Tj) and 
(f3"2"| . which then yield critical values fl s t(m) of f2, above 
which all the vortex solitons are azimuthally stable for 
charge to. These values can then be compared to those 
computed by the numerically exact profiles. 



C. Analytical Stability Predictions Using the 
Variational Profile 

As we showed in Sec. IIIH the variational ansatz yields 
a useful approximation to the true radial profiles. There- 
fore, the VA can be employed to calculate the constants 
in Eq. (|12p . and thereby derive analytical predictions for 
the azimuthal modulational stability. In Sec. IIVBI it was 
demonstrated that, for studying the stability of the vor- 
tices, only the values of -Kcrit and C3/C1 are required, 



K, 



'60* - (15/8) + 5%/3fF/ (AT) 



crit 



n* - G(o*) 



9i 



G(n*) 



(33) 



(34) 



since the other coefficients can be expressed in terms of 
these two (here, as before, f2 = G(Q*)). 

The VA can predict the critical value of O* above which 
the vortices are modulationally stable. Inserting the con- 
stants defined in Eq. (TT2"j) into the stability criteria (f3~Tj) 
and l|32p . we arrive at the following expressions which de- 
termine the critical frequency: 



K. 



crit 



: 



= 0; 



(35) 
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predicts that azimuthally stable vortices exist for all val- 
ues of m, at 17 > fijf. We also solved Eq. (p6|) for 17 at 
various values of to. The results are displayed in Fig. [5] 
It is seen that the critical frequency S7™(m) increases with 
the increase of to and eventually converges to fl™. Thus, 
the stability window in f2 is larger for lower charges. In 
Rcf . [27j , it was also concluded that azimuthally stable vor- 
tices exist for all to, and that lower charges indeed have a 
larger stability window. However, according to Ref. [27) . 
fist = £^™ x , i-e., the stability window shrinks towards as 
to — y co. The estimate obtained in that work shows that 
the window shrinks as 1/m 2 . This conclusion precludes ex- 
perimental creation of higher-order stable vortices, as the 
respective stability interval would be too small, and the 
radius of the vortex too large. 

According to the VA predictions, fi s t does increase with 
to (and, as shown in Fig. [5J the difference fi s t — f2 s t(w) 
decays proportionally to 1/to 2 , which resembles the pre- 
diction of Ref. [27[). On the other hand, since the varia- 
tional^ predicted value of il s t is smaller than 3/16, there 
remains a stability window at all to, which does not vanish 
at to — > co. 

As shown in Table [TT1 the variational predictions are very 
close to those using the numerically-exact profiles. We then 
note that the numerical results reported in Ref. [27j are 
more accurate than our predictions, the stability window 
indeed shrinking to zero at high values of to. 



AZIMUTHAL MODULATIONAL STABILITY: 
NUMERICAL RESULTS 



FIG. 8: (Color online) Top: The critical value of frequency 
f2™(m) versus topological charge m, according to the prediction 
of the variational approach. Bottom: The rate of the conver- 
gence between f2™(m) and f2™ at increasing m. The plotted 
curve, const x m~ 2 , which starts from our m = 1 computed 
value, demonstrates that the convergence rate is proportional 
to 1/m 2 . 



A', 



'6ft* 



crit 



- = o, 

TO 



(36) 



where T is defined as per Eq. ([13]) . We see that, at to — > co, 
these two criteria tend to coincide. At finite to, Eqs. (|35]l 
and (|3"6"|) give different critical values of il* , and hence of Q, 
too — one which depends on charge to, and the other one, 
independent of to, which represents an upper bound on 
-RT crit for all values of to. We denote the charge-dependent 
critical value as J7 s t(m), and the charge- independent one 

as n st . 

Solving Eq. (|35|) for f2* by means of a root finder and 
inserting the result into Eq. (|15[) yields 



(37) 



= 0.144320424. 



Since this value is smaller than fi^x = 3/16, this result 



In this section we report the numerical predictions 
and full simulation results for the azimuthal modulational 
(in)stability of the vortices. For the 2D simulations, we 
used a finite-difference scheme with a central difference in 
space and fourth-order Runge-Kutta in time [l6| . We used 
both polar-coordinate and Cartesian grids. The polar grid 
makes computing the growth of individual perturbation 
modes easier, therefore we used it to test our predictions 
for azimuthally unstable vortices. However, the polar grid 
forces one to use smaller time steps in the finite-difference 
scheme than is needed for the equivalent Cartesian grid. 
Therefore, for testing the critical frequency D, s t(m), which 
requires running multiple long-time simulations, we use the 
Cartesian grid. We used second-order differencing for the 
polar-coordinate simulations, and (due to the long simula- 
tion time needed) a fourth-order differencing for the Carte- 
sian simulations. 



A. Unstable Vortices 

In Fig. [3] we show the results for the (in)stability of vor- 
tices with charges to = 1,...,5 and ft <G [0.03,0.14]. In 
general, we see a good agreement between the numerically 
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FIG. 9: (Color online) Numerical predictions and numerical results for the growth rates of the fastest growing unstable perturbation 
mode around vortex solitons with charges m = 1, ...,5 (left to right, top to bottom). In the interval of Q £ [0.03,0.14], with the 
step of AQ = 0.01, the growth-rate predictions (shown by circles) were obtained by computing the integrals in the underlying 
Lagrangian, using the numerically exact profiles obtained through the Gauss-Newton optimization. We then ran full simulations of 
the vortex and recorded the average instability growth rate (shown by squares). The radial, angular, and temporal variables were 
discretized with steps Ar = 1, A6 = (2n)/(20 max[m, K max , 2]), and At = 0.001, respectively. Overall, the predictions match the 
numerical results very well. 



measured growth rates and the predictions for m > 2; how- 
ever, for values of which get closer to 3/16, the accuracy 
of the predictions becomes low in each case, implying that 
our predictions for f2 st (m) are not precise. To further test 
this, we ran long-time simulations of randomly perturbed 



vortices. 

For m = 1,2,3, we were easily able to identify the tran- 
sition from unstable to stable vortices, but for m > 3 we 
found it very difficult to detect a stable solution, because of 
a snake-like instability which breaks the vortex into asym- 
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FIG. 10: (Color online) Examples of the evolution of a vortex 
with charge m — 4 and fi = 0.18, perturbed by a random 
perturbation of size e = 0.05. Pictured is the squared absolute 
value of the wave function. We have set Ar = 2, r ma x = 120 
and At = 0.6. The initial vortex shape becomes deformed and 
then breaks up irregularly. 



metric irregular fragments. An example of this is shown 
in Fig. [TUJ Obviously, this snake-like instability is not cap- 
tured by our study of the azimuthal modulational instabil- 
ity. It is very plausible that it accounts for the discrepancy 
between the analytic and numerical growth rates observed 
in Fig. M 



B. Stable Vortices 

To identify the values of il s t (m) corresponding to the sta- 
bility border, we simulated the evolution of vortices with 

taken near the predicted values of 51 st (m), perturbed 
by a small uniformly distributed random noise in the az- 
imuthal direction. Our aim was to find a value Oi of Q that 
results in an unstable vortex soliton (and to observe the ac- 
tual breakup), and then show that the vortex solution for 
SI2 = ^1 + 0.001 is stable, by simulating its evolution for 
an extremely long time, compared to the time necessary to 
reveal the full breakup in the former case. We did this for 
various charges m. In Fig. II H we show the initial and final 
states produced by this analysis for m = 1. 

Since the snaking effect hinders our ability to simulate 
the evolution of the vortices for extended time periods at 
m > 3, we were not able to make precise stability predic- 
tions in this case. However, as this effect is dynamically 
distinct from the azimuthal breakup, we can still give an 
approximate estimate of the critical frequency for higher 
charges. 

In Table HI] we display the predictions of f2 s t (m) and 
their numerically found counterparts for m = 1, 5. The 
column labeled "semi-numerical" are the predictions made 
using the numerically exact radial profiles, while the col- 
umn labeled "VA" gives the predictions made using the 
variationally derived profiles. For comparison, in the same 
Table we also display findings reported in earlier works. 
We see that most of the predictions for m = 1 are close to 
the numerical result. For m > 1, the results reported in 
Ref. " 



27j are most accurate ones. 
It is observed that the semi-numerical and variational re- 
sults are very close, the latter ones being, in fact, slightly 
more accurate in comparison with the numerical findings. 
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FIG. 11: (Color online) An example of the numerical identi- 
fication of the stability border, Q s t(m), for m = 1. In this 
case, the evolution of the vortex with fli = 0.146 (top left) was 
simulated with a random-noise perturbation in the azimuthal 
direction, until it broke up into fragments due to the azimuthal 
instability (top right). We then carried out long simulations 
of the evolution of a vortex with Q2 = fii + 0.001 = 0.147 
(bottom left), to make it sure that the vortex is stable — in 
this example, up to t = 50,000 (bottom right). The drift of 
the position of the vortical pivot is due to the momentum im- 
parted to it by the random perturbations. We thus conclude 
that Q s t(m = 1) 6 [0.146,0.147]. In this example, we set the 
grid spacing to be Aa: = Ay = 1, time step to At = 0.2, and 
the perturbation amplitude to e = 0.05. 



m 


NUM 


scmi-NUM 


VA 


Ref. [27] 


Ref. [23] 


Ref. [31] 


1 


0.147 


0.1399 


0.1403 


0.1487 


0.16 


0.145 


2 


0.162 


0.1433 


0.1434 


0.1619 


0.17 


NC 


3 


0.171 


0.1437 


0.1439 


0.1700 


NC 


NC 


4 


(0.178) 


0.1437 


0.1441 


0.1769 


NC 


NC 


5 


(0.18) 


0.1437 


0.1442 


0.1806 


NC 


NC 



TABLE II: Comparison of the predictions and numerical results 
for the stability border, tt s t(m). The columns labeled "semi- 
NUM" and "VA" represent the predictions utilizing the numer- 
ically exact and the variationally derived radial profiles respec- 
tively. The column labeled "NUM" are the results obtained 
from the full 2D simulations. Numerical findings in parenthe- 
ses are those that were hard to fix due to the emergence of the 
snake-like instability (not comprised by our stability analysis), 
and may therefore be less accurate than the others. We also 
show predictions reported in earlier works. When no value has 
been computed or reported, we label the entry as "NC". The 
predictions made in Ref. [27J are the most accurate ones, in 
comparison to our simulations. 



This implies that the discrepancies in the variationally de- 
rived predictions are not due to any inherent deficiency in 
the VA. Rather, the discrepancies of both sets of predic- 
tions are likely due to the approximation of the vortex soli- 
ton as a separable entity composed of radial and azimuthal 
parts (as discussed in Sec. IIV[) . 
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FIG. 12: (Color online) Collisions between stable vortices, all with frequency Q — 0.17. First column: an elastic collision between 
a moving vortex of charge m — 2, with a slow initial velocity of 0.03, and a stationary vortex of charge m — 1. The vortices have 
no relative phase difference between them, which corresponds to their colliding sides being out-of-phase, hence the interaction is 
repulsive. Accordingly, the vortices undergo an elastic collision. Second column: the same as in the first column, but with a n phase 
difference between the colliding vortices, which implies that their colliding sides are in phase. This, in turn, gives rise to a mutual 
attraction, which causes the vortices to merge and eventually break up into fragments. Third column: same as in the first column, 
but for a larger collisional velocity of 0.3. The mutual repulsion between the in-phase vortices is not enough to counterbalance 
the high collisional momentum, and the vortices merge and break up into fragments. Fourth column: a charge m = 1 vortex with 
velocity 0.03 undergoing an elastic collision with an m = 2 vortex, in the case of a ir phase difference between them, which in 
turn collides with another 7r-phase-shifted vortex of charge m — 1. This phase arrangement corresponds to vortices colliding with 
adjacent sides that are mutually 7r-out-of-phase, and thus repel each other. These results attest to the robustness of the vortices. 



VI. COLLISIONS AND SCATTERING OF 
STABLE VORTICES 



In this section we study collisions and scattering of stable 
vortices in the CQNLS. In Rcf. [3l[, such collisions of vor- 



tices of unit charge were considered in a brief form. It was 
shown that both elastic and destructive collisions could be 
observed depending on the phase difference, A<f>, between 
the colliding vortices, as well as on their relative velocity. 
Here we expand this study in three directions: i) we in- 
clude vortices of charge m = 2; ii) we consider the critical 
velocity for elastic collisions as a function of A</>; iii) we 
study the scattering of vortices colliding at various values 
of impact parameters, by measuring the scattering angles 
and speeds. 

Summarizing results of numerous simulations, we have 
found that, at sufficiently small collisional velocities [39|. 
the elastic or destructive character of the collision is solely 
determined by the phase difference at the point of contact. 
This fact seems to be independent of the charge of the vor- 
tices involved in the collision. In Fig. [12] we display four 
different cases that illustrate the main features of head-on 
collisions between vortices with different charges. In the 
first two columns of Fig. Q2] we show the collision of vor- 
tices carrying charges m = 2 and m = 1, with A(f> = 
and A(j> = it, respectively. It can be seen that, since the 
vortex with m = 2 has opposite phase on the collision side, 
with respect to the vortex with m = 1, the interaction is 
repulsive since it locally emulates the interaction of two 



out-of-phase fundamental bright solitons 0, [2l|. There- 
fore, if the velocity is small enough, the mutual repulsion 
determines the result of the collision. On the other hand, 
when the vortex with charge m = 2 is phase-shifted by 
7r (or similarly, if the charge m = 2 vortex were to be 
placed on the opposite side of the charge m = 1 vortex), 
the adjacent sides of the colliding vortices are in phase. 
Therefore, the interaction is similar to that between two 
in-phase bright solitons, which is attractive [l], [2l[. This 
results in the merger of the two vortices and their eventual 
breakup into fragments. 

Clearly, even when the repulsive force acts between vor- 
tices, but the collision velocity is sufficiently large, the vor- 
tices have enough momentum for causing them to merge 
and eventually break up as well. This case is depicted in 
the third column of Fig. Q~2J 

Slowly moving vortices with opposite phases at the point 
of contact are quite robust against the collisions, as illus- 
trated in the fourth column of Fig. [T^l In this panel we 
depict a "billiard" -type example, in which one vortex col- 
lides with another, which in turn collides with a third one. 
All collisions in this case are elastic because the relative 
phases have been chosen such that the colliding sides are 
always out-of-phase, providing for the necessary repulsion. 
It is worthy to note that, although these collisions at low 
velocities seem to be elastic, the stable vortices clearly have 
internal breathing modes, that are excited during the col- 
lisions. Therefore, a small fraction of the collisional energy 
is transferred to these internal modes, preventing the colli- 
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FIG. 13: (Color online) Critical velocity w C rit for elastic collisions between two vortices as a function of their phase difference A<f>. 
The vortex charges of the colliding vortices correspond to: (a) (mi, HJ2) = (1, 1), (b) (mi, 7712) = (1, 2), and (c) (mi, 7712) = (2, 2). 
All vortices are taken with intrinsic frequency f2 = 0.17. 
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FIG. 14: (Color online) Scattering between an incoming vortex with charge mi and velocity v\ nl , and a quiescent vortex of charge 
7712. The trajectory for the incoming (stationary) vortex is depicted by small circles (squares). All panels depict different trajectories 
for different values of impact parameter d, as indicated [d = 5, 15, 25, 65 from top to bottom (bottom to top) for the incoming 
(stationary) vortex]. Panels (a) and (b) correspond, respectively, to two unitary-charge vortices (mi = m2 = 1), with A(f> = 
and collision velocities v\ nl = 0.03 and vf 11 = 0.06. Panels (c) and (d) correspond, respectively, to a mi = 1 vortex scattered by a 
7712 = 2 one and vice versa, for the collision velocity v\ nl = 0.03 and A<f> = 7r/2. Panel (e) corresponds to a mi = 2 vortex scattered 
by a ?7i2 = 2 one, for v\ nl = 0.02 and A(f> — n/2. All vortices have intrinsic frequency Q = 0.17. 



sions from being completely elastic. A more in-depth anal- 
ysis of the excitation of these internal modes and of the 
degree of the inelasticity of the collisions fall outside of the 
scope of the present work. Nonetheless, in order to quantify 
the degree of the elasticity of the collisions, we have mea- 
sured the critical velocity v cr it for the vortices to elastically 
bounce off each other: at v < 7> cr it the vortices bounce elas- 
tically, while at v > v CT n they merge and, most often, get 
destroyed. The crucial parameter that controls v cr it, for a 
particular combination of the charges of colliding vortices, 
is their phase difference A</>. Different values of A(j> allow to 
have different relative phases at the colliding sides (as dis- 



cussed above), and thus (partially) repel or attract. This 
effect is illustrated by Fig. 1131 where we depict w cr i t versus 
A(f) for vortex pairs with charges (a) (mi, 7712) = (1, 1), (b) 
(7711,7772) = (1,2), and (c) (7711,7772) = (2,2). From panel 
(a) it is clear that, for charges mi = TO2 = 1, the most "sta- 
ble" (in terms of observing elastic collisions at larger colli- 
sional velocities) collision occurs around A</> = 0. This cor- 
responds exactly to the most stable case mentioned above, 
when the collisional sides are out of phase. In panel (b) we 
depict v cr n for (mi, 7712) = (1,2). It is seen in this panel 
that v cr it is 7r-periodic, as a function of A</>, due to the 
angular symmetry of the 1112 = 2 vortex. For this vortex- 
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FIG. 15: (Color online) Scattering angles for the collision of two vortices corresponding to the panels depicted in Fig. 1141 Depicted 
by (blue) circles is the negative of the scattering angle —9i of the incoming vortex, while (green) squares stand for the scattering 
angle 62 of the initially quiescent vortex. Triangles (in red) display the difference between these two angles, 62 — 61. All angles are 
given in degrees. 



charge combination, the most stable collision occurs around 
A<fi = 7r/2 (or, due to the symmetry, A0 = 3w/2), which 
again corresponds to the adjacent out-of-phase sides. Fi- 
nally, in panel (c) we present the results for two vortices 
with mi = m2 = 2. The behavior is very similar to that 
in panel (b): the same periodicity, shape and value of the 
optimal phase difference. However, since the vortices with 
m = 2 vortices are less stable than those with m = 1, 
the critical velocity u cr ;t is lower for the (mi,m,2) = (2,2) 
charge combination than for the case of (mi, 7712) = (1, 2). 

In order to fully characterize the collisions between vor- 
tices with different charges, we studied the scattering of a 
vortex of charge mi impinging at velocity vf 11 upon another 
vortex of charge m2 at rest (v 2 m = 0), with impact parame- 
ter d. The impact parameter d is defined as the perpendic- 
ular distance between the trajectory of the incoming vortex 
and the initially quiescent one. In Fig.Q3]we show several 
orbits for the two vortices with different impact parame- 
ters, vortex-charge combinations, and incoming velocities. 
It is seen, in all the panels, that all the trajectories are 
well defined for the chosen parameters, and clearly feature 
elastic scattering. The key to achieve elastic scattering is 
to choose, for each charge combination, the optimal phase 
difference, using Fig. [T3] [i.e., A</> = for mi = m 2 = 1 
and A<f) = tt/2 for (mi,m 2 ) = (1,2), (2, 1), or (2,2)], and 
the initial velocity lower than the corresponding i> cr it. For 
example, we chose in panel (a) of Fig. [TJ] A0 = and 
4 ni = 0.03 < v C rit(7T = 0) = 0.067. In fact, when one 
uses a velocity closer to the critical one the vortices start 
experiencing strong deformations of the quadrupole and oc- 



tupole types when they collide nearly head-on. This effect 
is clearly visible for the trajectory with d = 5 (also d = 10 
and d = 20) in panel (b) of Fig. M where v 2 ni = 0.06 is 
close to the critical value, w cr i t = 0.067, and the trajectory 
features undulations due to the internal deformation of the 
vortex. Even larger velocities, i.e., above i> cr it, lead to elas- 
tic scattering for larger values of the impact parameter d, 
and to annihilation of the colliding vortices at smaller d 
(not shown here). It is interesting to note that, as for hard 
spheres with different masses, the scattering of vortices can 
also produce very different scattering angles, depending on 
the mass combination of the two colliding objects. This 
effect is made evident by the comparison of the scattering 
between vortices with mi = 1 and m-2 = 2 in panel (c), and 
the scattering between mi = 2 and m2 = 1 in panel (d) 
(all the other parameters coincide in both panels). When 
a "lighter" vortex with mi = 1 collides with the "heavier" 
one, with mi = 2, for d < 15, the former vortex bounces 
back after the collisions. In the reverse situation, with the 
"heavier" TO2 = 2 vortex impinging upon the "lighter" one 
with mi = 1, the "heavier" vortex only loses a portion of 
its momentum and continues to move in a straighter trajec- 
tory, while the "lighter" vortex is pushed away at a larger 
velocity. 

To further characterize the scattering between vortices 
we analyzed the scattering angle as a function of the im- 
pact parameter. In our scattering experiments, we initially 
define horizontally moving and quiescent vortices. After 
the collision, the trajectory of the moving vortex is de- 
viated by scattering angle 9\ (in our case negative since 
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FIG. 16: (Color online) Final (post-collision) velocities for the scattering of two vortices corresponding to the panels depicted in 
Fig. HU Depicted in (blue) circles and (green) squares are, respectively, the final velocities, vi and V2, of the initially moving and 
quiescent vortices. (Red) triangles depict the sum of these two velocities, Vi + i>2- The thin horizontal black line designates the 
initial velocity of the moving vortex. 



we consider the moving vortex placed below the quiescent 
one), and the initially quiescent vortex is scattered at angle 
6*2 (in our case, positive) . In Fig. [15] we depict the scatter- 
ing angles corresponding to the scattering experiments dis- 
played in Fig. [Mj We depict by (blue) circles the negative 
of the scattering angle, —B\, of the incoming vortex, and 
by (green) squares the scattering angle 6* 2 . In all the cases, 
we also depict the difference between the scattering angles, 
A0 = 02 — 01, by (red) triangles. In the case of the elastic 
collision between two hard spheres of equal masses it is well 
known that A0 = 90°. Figure PT~51 reveals that the collisions 
between vortices with equal charges [mi = 777,2 = 1 in pan- 
els (a) and (b), and mi = 7712 = 2 in panel (e)] produce 
A0 very close to 90° for all values of the impact parameter. 
In contrast, as it is the case for hard spheres with differ- 
ent masses, our simulations demonstrate that A0 > 90° 
for mi < ?772 [(mi,m 2 ) = (1,2) in panel (c)], and that 
A0 < 90° for mi > m 2 [(mi, m 2 ) = (2, 1) in panel (d)]. 

Another measure of the scattering can be given in terms 
of the initial and final velocities. In Fig. [THl we present 
the final velocities corresponding to the scattering events 
shown in Fig. 1141 The figure shows, by (blue) circles and 
(green) squares, the post-collision velocities, v\ and V2, re- 
spectively, of the initially moving and quiescent vortices. 
We also depict, using (red) triangles, the sum of the final 
velocities, v\ + t> 2 , and the thin horizontal black line repre- 
sents the initial velocity of the moving vortex, v™. It can 
be concluded that v\ + V2 < vf 11 for relatively small values 
of the impact parameter, d, and v\ + tj 2 > vf 11 for larger 



values of d. 



The above results suggest that the interactions between 
vortices may be either elastic or destructive, depending on 
the parameters. In order to simulate the interplay of these 
interactions in a multi- vortex gas, we display in Fig. [T7] a 
time series depicting the time evolution of a random vortex 
gas. The first two rows depict the local intensity plots at 
the indicated moments of time, while the bottom two rows 
depict the corresponding field phases. We use a square do- 
main of size 120 x 120 with periodic boundary conditions, 
seeding 15 vortices of charge m = 1 at random positions 
(such that initially no two vortices were closer than twice 
the size of a single vortex), and with random phases. As 
seen from the figure, due to the random phase differences 
and collision angles, some vortex pairs experience elastic 
collisions, while others tend to merge and destroy each 
other. Since the destruction of vortices is an irreversible 
process, all the vortices get eventually destroyed, leading 
to a seemingly disordered pattern of interacting intensity 
patches with an approximately constant phase inside each 
one (see the correspondence between field patches in the 
top rows and the phase patches in the bottom rows, at the 
late stage of the evolution). The dynamics resembles the 
grain coarseningin typical chaotic spatiotemporal systems, 
see e.g., review J3(. 
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FIG. 17: Interaction of a random gas of vortices of charge in = 1 in a periodic box. The initial separation between any two 
neighboring vortices exceeds twice the width of a single vortex. The initial phases are set randomly. The first two rows depict the 
evolution of the field intensity at the times indicated, while the bottom two rows depict the respective phase of the field. Note the 
coarsening effect due to vortex merger and destruction. 



VII. CONCLUSIONS 

In this work, we have revisited the issues of the existence, 
stability, and interactions of solitary vortices in the two- 
dimensional CQNLS (cubic-quintic nonlinear Schrodingcr) 
equation, using both numerical and analytical methods. 
The latter one is based on the variational approach, that 
has been shown to be increasingly more accurate as the vor- 
tex' topological charge increases. We have also developed 
the analysis of the azimuthal modulational (in) stability of 
vortices, using the approximation originally proposed in 
other contexts in Refs. Q and 0, which postulates the 
separation between the frozen radial profile of the vortex 
soliton and free evolution of azimuthal perturbations. This 
approach leads to an effectively one-dimensional equation 
for the azimuthal dynamics. Examining the stability of the 
perturbed annulus within the framework of the latter equa- 
tion, we were able to predict the stability of the vortices 
against the breakup. We then ran full 2D simulations to 
verify the semi-analytical predictions. 

For azimuthally unstable vortices, our predictions of the 
largest growth rate are fairly accurate over a wide range 
of the vortex' intrinsic frequency (especially for the higher- 
order vortices, with topological charge m > 2). The so de- 
veloped semi-analytical technique may be useful for other 
applications. The analytical approach is less accurate in 
predicting the stability border (intrinsic critical frequency 



of the vortex) for m > 1. This discrepancy is due to de- 
viation from the assumption admitting the separation of 
the perturbed evolution into radial and azimuthal parts, 
while the VA proper, which was used to predict the radial 
shape of the stationary vortex solitons turns out to yield 
accurate results. Comparing our numerical results for the 
critical frequency with previously published ones, we have 
concluded that the most accurate findings were reported 
in Ref. 27 1 . We have also found that higher-order vortices 



may be subject to a snaking radial instability. 

The semi-analytical approach developed in this work 
may be quite useful for other 2D models — at least, for the 
description of the most stable vortex solitons with lowest 
values of the topological charge. 

We have also investigated, in detail, collisions between 
stable vortex solitons in the two-dimensional CQNLS equa- 
tion. The outcome of the collision crucially depends on the 
phase difference and relative velocity of the vortices. We 
have produced the results showing the critical velocity for 
elastic collisions as a function of the relative phase for all 
combinations of the vortices with m = 1 and 2. We also 
studied the dependence on the scattering angles and post- 
collision velocities on the impact parameter of the collision. 

One relevant direction for the extension of the present 
analysis is to find accurate conditions for the stabiliza- 
tion of the nearly isotropic vortex solitons by an external 
spatially periodic (lattice) potential, see Ref. [33| and ref- 
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erences therein. A challenging problem is to extend the 
semi-analytical considerations to three-dimensional soli- 
tons with embedded vorticity, in the model with the same 
cubic-quintic nonlinearity. In fact, it is not known if such 
3D vortex solitons with higher vorticity (m > 2) may be 
stable [HHl. 
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